Monocyte-derived transcriptomes explain the ineffectiveness of abatacept in rheumatoid arthritis

Background The biological mechanisms underlying the differential response to abatacept in patients with rheumatoid arthritis (RA) are unknown. Here, we aimed to identify cellular, transcriptomic, and proteomic features that predict resistance to abatacept in patients with RA. Methods Blood samples were collected from 22 RA patients treated with abatacept at baseline and after 3 months of treatment. Response to treatment was defined by the European League Against Rheumatism (EULAR) response criteria at 3 months, and seven patients were classified as responders and the others as non-responders. We quantified gene expression levels by RNA sequencing, 67 plasma protein levels, and the expression of surface molecules (CD3, 19, and 56) by flow cytometry. In addition, three gene expression data sets, comprising a total of 27 responders and 50 non-responders, were used to replicate the results. Results Among the clinical characteristics, the number of monocytes was significantly higher in the non-responders before treatment. Cell type enrichment analysis showed that differentially expressed genes (DEGs) between responders and non-responders were enriched in monocytes. Gene set enrichment analysis, together with single-cell analysis and deconvolution analysis, identified that Toll-like receptor 5 (TLR5) and interleukin-17 receptor A (IL17RA) pathway in monocytes was upregulated in non-responders. Hepatocyte growth factor (HGF) correlated with this signature showed higher concentrations in non-responders before treatment. The DEGs in the replication set were also enriched for the genes expressed in monocytes, not for the TLR5 and IL17RA pathway but for the oxidative phosphorylation (OXPHOS) pathway. Conclusions Monocyte-derived transcriptomic features before treatment underlie the differences in abatacept efficacy in patients with RA. The pathway activated in monocytes was the TLR5 and IL17RA-HGF signature in the current study, while it was the OXPHOS pathway in the replication set. Elevated levels of HGF before treatment may serve as a potential biomarker for predicting poor responses to abatacept. These findings provide insights into the biological mechanisms of abatacept resistance, contributing valuable evidence for stratifying patients with RA. Supplementary Information The online version contains supplementary material available at 10.1186/s13075-023-03236-y.


Background
Rheumatoid arthritis (RA) is a chronic autoimmune and inflammatory disease that leads to joint destruction.Abatacept (CTLA4Ig) binds to CD80 and CD86 on antigen-presenting cells, thereby blocking the costimulatory interaction with CD28 on T cells and inhibiting its activation.This biological treatment is often considered when Phase I interventions, including the usage of methotrexate (MTX), do not sufficiently control disease activity of RA [1].Most of the patients receive beneficial effects from abatacept; however, some exhibit poor responses to the therapy [2,3], prompting extensive studies into factors influencing its outcome.
Numerous studies have investigated potential contributors to non-responsiveness.As the factors for nonresponders, low anti-citrullinated peptide (ACPA) titer [4,5], low rheumatoid factor (RF) titer [5,6], lower number of CD19 positive cells [7], higher number of CD28 negative T cells [8], higher interferon signature [9], and shared-epitope (SE)-negative [10] have been reported.However, no biomarkers have been considered superior to ACPA/RF that are clinically useful and biologically interpretable.Consequently, there is an unmet need for further exploration concerning this critical issue.
Here, we aimed to identify cellular, transcriptomic, and proteomic features that predict responses or resistance to abatacept in a Japanese RA cohort as well as using public gene-expression data.

Study design
We enrolled 22 RA patients who received abatacept therapy at the Center for Rheumatic Disease at Kyoto University Hospital.RA was diagnosed according to either the 1987 revised American College of Rheumatology classification criteria [11] or the 2010 America College of Rheumatology/European League Against Rheumatism (ACR/ EULAR) criteria [12].The overview of this study is shown in Fig. 1.Peripheral blood was collected from patients before and approximately 3 months after initiating abatacept treatment using heparin and ethylenediaminetetraacetic acid (EDTA) blood collection tubes.

Clinical evaluation and response measure
Disease activity was evaluated using the Disease Activity Score 28-Erythrocyte sedimentation rate (DAS28-ESR) at every clinic visit.Clinical characteristics, including age, sex, body mass index (BMI), smoking history, duration of RA, the titers of RF, anti-cyclic citrullinated peptide (CCP) antibodies, anti-nuclear antibodies (ANA), ESR, C-reactive protein (CRP) and treatment profiles [the use of MTX and prednisolone (PSL)] before the initiation of treatment, and white blood cell count, blood differential count before and 3 months after the initiation of treatment, were obtained from medical records.The Steinblocker stage was assessed by rheumatologists.We classified patients using the European League Against Rheumatism (EULAR) response criteria [13] at 3 months after treatment initiation.Moderate and no responders were classified as non-responders, and the others as good responders.

RNA sequencing, transcriptome analysis
RNA sequencing and quantification of gene expression were performed as described previously [14].In short, RNA was extracted from peripheral blood mononuclear cells (PBMCs), and sequencing was conducted by HiSeq 2500 in the 150-bp paired-end mode.We trimmed reads by Cutadapt (ver1.1),aligned to the GRCh37 reference genome by STAR [15] (ver.2.7.3a), counted gene expression by RSEM [16] (ver.1.3.1),performed normalization using size factor implemented in DeSeq2 [17], and converted to count per million (CPM).The Wald test performed a differential gene expression analysis using DESeq2.Gene set enrichment analysis [18] was performed for 32,880 gene modules in MSigDB (ver 7.5) [19] with 10,000 permutations.We focused on the results with positive enrichment scores.As for the analysis of the replication data set, gene set enrichment analysis was performed by Metascape [20] with an input of 166 genes contributing to monocyte enrichment.Cell-type enrichment analysis was performed by WebCSEA [21] with the default setting.In both the current study and the replication dataset, the top 500 genes with the lowest P values were used as input since this platform accepts no more than 500 genes as input.This method outputs enrichment statistics for each cell type, along with a list of genes contributing to the enrichment of each cell type.The HLA-DRB1 allele was typed using HLA-HD (ver.1.5.0)[22] and the shared epitope (SE) allele was classified as a previous report [23].

Flow cytometry, and protein measurements
Flow cytometry and protein measurements were performed as described previously [14].In short, we assessed the surface molecule expression by BD Canto ™ II using the following antibodies obtained from BD Pharmingen: allophycocyanin-conjugated anti-CD56 (341025, NCAM16); antiCD3 (566683, OKT3), and V500-conjugated CD19 (561121, HIB19).We gated lymphocytes based on forward scatter (FSC) and side scatter (SSC) parameters and then calculated the percentage of each cell fraction in lymphocytes.The analysis was conducted using FlowJo software.The absolute number of each cell fraction in peripheral blood (× 10 6 /mL) was calculated using the percentage and an absolute number of lymphocytes measured by the hematology analyzer MEK-7300 (Nihon Kohden ® ).We measured 67 proteins in plasma using ProcartaPlex Human 15-plex, Procarta-Plex Human 49-plex, Human VCAM-1 Simplex, Human sICAM-1 Simplex, and Human sCGF Simplex with Bio-Plex 200 (BIO-RAD) according to the manufacturer's instructions.These assays were selected during the study's design phase to maximize the number of measurable proteins through a multiplex system.

Replication data set
We searched the gene expression data concerning abatacept response in gene expression omnibus (GEO) as much as possible in January 2023.We utilized three datasets; the breakdown is shown in Table S1.In each case, responders and non-responders were defined by the available information.As for GSE78068 [24], responders were defined by remission according to the Clinical Disease Activity Index (CDAI) at six months of therapy.As for GSE68215 [25], non-responders were defined as moderate / no response according to EULAR criteria at six months of therapy.As for GSE172188 [26], non-responders were defined as moderate / no response according to EULAR criteria at 16 weeks of therapy.We performed an association test using limma (ver.3.46.0)[27] in each cohort and meta-analyzed the results by MetaVolcanoR (ver.1.14).The scripts to analyze each data set and meta-analysis are provided in the GitHub repository: https:// github.com/ takes hiiwa saki/ abata cept.We defined the oxidative phosphorylation (OXPHOS) score as the mean expression of 200 genes included in "HALLMARK_OXIDATIVE_PHOSPHO-RYLATION" in MSigDB (ver.7.5) [19].

Single-cell analysis and deconvolution analysis
To elucidate the cellular origin of the detected module, we used the single-cell data of PBMCs [28] and synovium [29] from RA patients.The process of annotating cell types was documented in the provided GitHub repository: https:// github.com/ takes hiiwa saki/ abata cept.We performed a deconvolution analysis using CIBERSORTx [30] to estimate the expression of genes in monocytes.As a reference, we used the single-cell RNA expression data of PBMCs from RA patients [28].

Ethics approval and consent to participate
The present study was performed following the Helsinki Declaration and was approved by the Kyoto University Graduate School and Faculty of Medicine Ethics Committee (approval number: G0511).Written informed consent to participate in the present study and publish the results obtained was provided by all enrolled patients.

Comparison of clinical characteristics of responders and non-responders
Twenty-two RA patients were enrolled in this study (Fig. 1).Seven (31.8%) patients had a good response to abatacept, while 15 (68.2%) did not.The clinical characteristics of patients at the treatment initiation are summarized in Table 1.There was no difference in baseline clinical characteristics between responders and nonresponders (P > 0.05).However, when we investigated cellular phenotype, we found a higher number of monocytes (P = 0.02) and a lower number of CD3 + cells with suggestive significance (P = 0.09) in non-responders before treatment.After 3 months of treatment, there remains a tendency of a higher number of monocytes and a lower number of CD3 + cells in non-responders (Table S2) (P = 0.09, 0.03, respectively).

Relationship between gene expression and response
We performed a differential gene expression analysis between responders and non-responders before treatment to identify transcriptomic features that predict response to abatacept.We identified 952 differentially expressed genes (DEGs) associated with treatment responses (false discovery rate; FDR < 0.05) (Fig. 2A, Table S3).After 3 months, the number of DEGs decreased to 19 (Fig. 2B, Table S4).We performed a gene-set enrichment analysis [16] using the transcriptomic data before treatment to characterize transcriptomic differences.As a result, we found the strongest enrichment in "module 1" (P = 1.0 × 10 −4 , Enrichment score = 0.89) (Fig. 2C, Table S5).We interpreted this enrichment based on the current transcriptomic data and existing literature (Supplementary Data S1) and determined it was attributed to the upregulation of Toll-like receptor 5 (TLR5), myeloid differentiation primary response protein 88 (MYD88), and interleukin-17 receptor A (IL17RA) in non-responders.These three genes are significantly correlated with each other in 44 specimens (22 samples with two-time points) (Fig. 2D, Table S6) and exhibited higher expression levels in non-responders both before treatment as well as 3 months after treatment (Fig. 2E).
Next, to detect the cellular origin of the upregulated transcriptomic feature in non-responders, we performed a cell-type enrichment analysis.The results showed the strongest enrichment on monocytes among all the major immune cell types in PBMCs (P = 8.7 × 10 −212 , Fig. 2F).When we estimated the expression of TLR5, MYD88, and IL17RA in monocytes by deconvolution analysis, the non-responders had higher expression before treatment initiation as well as 3 months after treatment (Fig. 2G), MYD88 showing significant upregulation in nonresponders before treatment as well as 3 months after treatment (Table S7).

Cellular origin of TLR5, MYD88, and IL17RA
To further investigate the cellular origin of the TLR5, MYD88, and IL17RA, we analyzed single-cell data from PBMCs of RA patients.First, we found that TLR5, as a representative gene of TLR5 and MYD88, and IL17RA are mainly expressed by monocytes (Fig. 3A, B).Furthermore, we identified a subset of monocytes that co-express TLR5 and IL17RA (42 out of 1298 monocytes; 3.2%), with their expression showing a correlation (Fig. 3C).Hereafter, we will refer to these monocytes as "TLR5 + IL17RA + monocytes." We also validated the cellular origin of TLR5, MYD88, and IL17RA in RA synovium using single-cell data.Consistent with the findings in PBMCs, monocytes were the immune cell type with the highest expression of TLR5 and IL17RA (Fig. 3D, E).We also found TLR5 + IL17RA + monocytes, and their expression exhibited a strong correlation (Fig. 3F).The proportion of TLR5 + IL17RA + monocytes among monocytes was 14.5% (87 out of 602 cells), which was 5.1-fold higher than that of PBMCs (P = 8.2 × 10 −18 , Fisher's exact test).

Identifying proteins associated with the detected signature and treatment response
To identify the protein associated with this transcriptomic signature and treatment response, we calculated the associations between the mean expression of the three genes (TLR5, MYD88, and IL17RA) and the 67 proteins measured by multiplex immunoassay (see Methods) using 44 specimens (22 individuals × 2 time points).Among the 67 protein levels, C-X-C Motif Chemokine Ligand 10 (CXCL10), IL-8, hepatocyte growth factor (HGF), and IL-20 showed significant association with the expression of this module (Fig. 4A, Table S8).Among the four proteins, we found HGF was significantly higher in non-responders before treatment (P = 0.04, Fig. 4B, Table S9).We confirmed that the expression of the HGF was significantly correlated with those of TLR5, MYD88, and IL17RA using transcriptomic data derived from 44 specimens (Fig. 4C, Table S6).To confirm the potential to predict the response to abatacept by measuring the levels of HGF before treatment initiation, we performed a receiver-operating-characteristic (ROC) analysis.The area under the ROC curve (AUC) for HGF was 0.78, and the most accurate cut-off level was 176.21 pg/mL (Fig. 4D).These results suggest that the level of HGF protein is associated with the detected transcriptomic feature as well as a good predictor of response.
We then investigated the cellular origin of HGF using the same single-cell data.We found that monocytes primarily express HGF (Fig. 4E).We further compared the expression of HGF between TLR5 − /IL17RA − monocytes and TLR5 + IL17RA + monocytes (Fig. 4F).As a result, we observed an enrichment of HGF expression in TLR5 + IL17RA + monocytes (4 out of 42 cells vs. 52 out of 9469 cells; 17.4-fold P = 1.4 × 10 −4 , Fisher's exact test).We did the same analysis using the synovium single-cell data.Consistent with the result in  PBMCs, HGF was expressed by monocytes (Fig. 4G) and was enriched in TLR5 + IL17RA + monocytes compared to TLR5 − /IL17RA − monocytes (Fig. 4H) (23 out of 87 cells vs. 87 out of 515 cells; 2.4-fold enrichment, P = 1.1 × 10 −4 ).When we calculated the levels of HGF and the number of monocytes in the current data (22 individuals × 2 time points), we observed a positive correlation (Fig. S1), suggesting the number of monocytes in peripheral blood could be an alternative biomarker for HGF.Furthermore, based on several previous literatures indicating the association between HGF and bone damage [31,32], we analyzed the relationship between the Steinblocker stage and the HGF levels before treatment, but found no significant association (Fig. S2).

Replication of the current results
We assessed the obtained results across various datasets, compiling a gene expression dataset related to abatacept efficacy to maximize our sample size.Through meta-analysis of the three data sets (Table S1), we identified 839 genes exhibiting nominal significance (P < 0.05, Fig. 5A, Table S10).As for TLR5, MYD88, and IL17RA, they did not show any significant differences between responders and non-responders (P > 0.05).However, differentially expressed genes were notably enriched for those expressed in monocytes (P = 2.9 × 10 −90 ) (Fig. 5B).The set of 166 genes that contribute to monocyte enrichment (Table S11) demonstrated enrichment for genes involved in the aerobic respiration / OXPHOS pathway (Fig. 5C, Table S12).To quantify this difference, we calculated the mean expression of the OXPHOS pathway (designated as the "OXPHOS score") and compared it between responders and non-responders in each study.A significant increase in OXPHOS score in non-responder was observed in the GSE68215 dataset (Fig. 5D), while there was no statistical significance in all the other datasets (Fig. S3).Finally, using single-cell datasets, we investigated the cellular source of the OXPHOS signature within immune cells.We found that the highest expression of the OXPHOS signature was observed in monocytes, both within PBMCs (Fig. 5E) and within the synovium (Fig. 5F).

Discussion
We have identified that HGF-producing TLR5 + IL17RA + monocytes play critical roles in the ineffectiveness of abatacept in our study.Although this finding was not replicated in the replication set, we have demonstrated that distinct monocyte-derived transcriptome features The pathway activated in monocytes was the TLR5 and IL17RA to HGF pathway in the current study.According to the single-cell data derived from PBMCs and synovium from RA patients, there were TLR5 + IL17RA + monocytes, and HGF expression was enriched in that cell type.IL17A is known to modulate both expression of TLR5 and IL17RA [33,34].Expression of IL17A did not show a significance difference between responders and non-responders (Tables S3, S4).This lack of significance may be attributed to the limited statistical power arising from the low expression of IL17A (Tables S3, S4).Therefore, increasing RNA sequencing coverage could potentially capture the variations in IL17A expression and unravel the intricate relationship between its expression, IL17RA, and TLR5.Furthermore, we noted an association of HGF with TLR5, MYD88, and IL17RA signatures.HGF is associated with osteoclastogenesis in Collageninduced arthritis (CIA) mice [31] and radiographic damage in RA patients [32].Consequently, our findings gain significance by connecting abatacept treatment resistance to HGF's acknowledged pathological role in RA.
The OXPHOS is a metabolic pathway to produce adenosine triphosphate (ATP) in mitochondria and is involved in many immune cell functions [35] and dysregulation of this pathway is related to multiple autoimmune diseases, including RA [36].It is reported that increased levels of this signature in monocytes in RA [37].Our study introduces the intriguing possibility of utilizing this OXPHOS signature within monocytes as a tool to stratify RA patients.
This study has the limitation of a relatively small sample size in the current data.Especially, the re-evaluation of the cut-off level for HGF (Fig. 4D) is warranted in an independent and larger cohort.Nevertheless, our study yields valuable insights with direct clinical implications.Firstly, it is imperative to exercise caution when dealing with RA patients exhibiting elevated HGF levels, as this has been linked to poor responses not only to abatacept but also to TNF inhibitors [14] and the advancement of radiographic damage [34].Although we couldn't identify a significant association between HGF and bone damage in the current dataset (Fig. S2), further analysis is warranted to explore the potential effect of HGF on longterm bone damage.Secondly, our findings suggest an approach for patients whose monocytes exhibit activation of the TLR5 and IL17RA to the HGF pathway or the OXPHOS pathway.In such cases, exploring treatment options targeting innate immune cells' modulation could offer a promising alternative strategy.Interventions such as anti-IL-6 or anti-GM-CSF or JAK inhibitor [38,39] may be particularly effective in these scenarios, potentially offering enhanced clinical outcomes.

Conclusions
We have demonstrated that monocyte-derived transcriptomic features before treatment underlie the differences in abatacept efficacy.In the current study, the pathway activated in monocytes was the TLR5 and IL17RA to HGF pathway, while in the replication set, it was the OXPHOS pathway.The levels of HGF before treatment initiation may serve as a potential biomarker for predicting poor responses to abatacept.

Fig. 2 Fig. 3
Fig. 2 Transcriptomic analysis of the current study.A, B Volcano plots showing differentially expressed genes between responders and non-responders before treatment (A) and three months after treatment (B).C Results of gene set enrichment analysis before treatment.D Correlation between each gene's expressions.The number in each box represents the Pearson's correlation coefficient.E Expression was standardized across all samples, and Z scores were shown.Sample names beginning with "R" indicate responders, "NR" indicates non-responders; sample names ending with "0 M" indicate pre-treatment, and "_3M" indicates 3 months post-treatment.F Results of cell type enrichment analysis before treatment.T, T cell; NK, NK cell; B. B cell; DC. dendritic cell; Mono, monocyte.G Estimated gene expression in monocytes.As in E, Z scores are shown

Fig. 4
Fig. 4 Identification of HGF as an alternative biomarker of TLR5-MYD88 and IL17RA pathway.A Before treatment of responders (blue) and non-responders (orange), three months after treatment of responders (green) and non-responders (red).The Pearson's correlation coefficient (r), and the P-value calculated from linear regression are shown.B Comparison of HGF levels before treatment data.The P-value calculated by the Mann-Whitney U test is shown.C Correlation between each gene's expression.The number in each box represents Pearson's correlation coefficient.D ROC curves for no response to abatacept.E-H Expression of HGF in immune cells using single-cell data from peripheral blood (E, F) and synovium (G, H).CD4, CD4 + T cell; CD8.CD8 + T cell; NK.NK cell; B. B cell; Mono, monocyte

Fig. 5
Fig. 5 Transcriptomic analysis of the replication set.A Transcriptomic association study between responders vs. non-responders (meta-analysis of three datasets).Genes related to the OXPHOS pathway are color-coded(red: upregulated in non-responders, blue: downregulated in non-responders).B Results of cell type enrichment analysis.C Enrichment analysis of genes contributing to monocyte enrichment.D The OXPHOS scores in the GSE68215 data set.The P-value calculated by the Mann-Whitney U test is shown.E, F The OXPHOS scores in immune cells in peripheral blood (E) and in synovium tissue (F) in patients with RA.CD4, CD4 + T cell; CD8, CD8 + T cell; NK, NK cell; B, B cell; Mono, monocyte

Table 1
Clinical and cellular characteristics of responders and non-respondersData were described as medians (interquartile range (IQR)) for continuous variables and numbers (percentages) for categorical variables OR Odds ratio, ND No data, Inf Infinity, BMI Body mass index, CRP C-reactive protein, ESR Erythrocyte sedimentation rate, RF Rheumatoid factor, CCP Cyclic citrullinated peptide, PSL Prednisolone, MTX Methotrexate, b/tsDMARDs biologic/targeted synthetic disease-modifying antirheumatic drugs, including etanercept, tocilizumab, infliximab, adalimumab, golimumab, cDMARDs conventional DMARDs, including methotrexate, bucillamine, salazosulfapyridine, tacrolimus, cyclosporine, mizoribine, leflunomide, and iguratimod a The Mann-Whitney U test was used for continuous variables and Fisher's exact test for categorical variables, †Calculated from all samples regardless of the corresponding drug usage